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Abstract 

We investigate modified steepest descent methods coupled with a 
loping Kaczmarz strategy for obtaining stable solutions of nonlinear 
systems of ill-posed operator equations. We show that the proposed 
method is a convergent regularization method. Numerical tests are pre- 
sented for a linear problem related to photoacoustic tomography and 
a non-linear problem related to the testing of semiconductor devices. 
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1 Introduction 

In this paper we propose a new method for obtaining regularized approxi- 
mations of systems of nonlinear ill-posed operator equations. 

The inverse problem we are interested in consists of determining an un- 
known physical quantity x E X from the set of data (yo, . . . ,yj^-i) G , 
where X, Y are Hilbert spaces and > 1. In practical situations, we do 
not know the data exactly. Instead, we have only approximate measured 
data yf & Y satisfying 

\\y! - ViW < Si , i = o,...,N-i, (1) 

with Si > (noise level). We use the notation 6 := {Sq, . . . ,Sn-i)- The finite 
set of data above is obtained by indirect measurements of the parameter, 
this process being described by the model 

F,{x)=yi, i = 0,...,N-l, (2) 
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where Fj : C X — > y, and Di are the corresponding domains of definition. 

Standard methods for the solution of system (2) are based in the use 
of Iterative type regularization methods [1, 7, 13, 16, 19] or Tikhonov type 
regularization methods [7, 23, 30, 32, 33] after rewriting (2) as a single 
equation F(x) = y, where 

N-l 

F:=(Fo,...,F;v-i) : f] ^ (3) 

and y := {yo, . . . ,yN-i). However these methods become inefficient if N 
is large or the evaluations of Fj(x) and F'-{x)* are expensive. In such a 
situation, Kaczmarz type methods [6, 15, 22, 26] which cyclically consider 
each equation in (2) separately are much faster [24] and are often the method 
of choice in practice. 

For recent analysis of Kaczmarz type methods for systems of ill-posed 
equations, we refer the reader to [4, 10, 11, 17]. 

The starting point of our approach is the steepest descent method [7, 29] 
for solving ill-posed problems. Motivated by the ideas in [10, 11], we propose 
in this article a loping Steepest-Descent-Kaczmarz method (l-SDK method) 
for the solution of (2). This iterative method is defined by 

^i+i = ^k- ^kOkSk , (4) 

where 



Sfc:=Fjfc](xf)*(F[fc](xf)-4]), (5) 




Here amin > 0, r G [2,cxd) are appropriate chosen numbers (see (13), (14) 
below), [k] := {k mod N) G {0, . . . , — 1}, and Xq = xq X is an ini- 
tial guess, possibly incorporating some a priori knowledge about the exact 
solution. The function <I>i.ci '■ (0, oo) — > (0, oo) defines a sequence of relax- 
ation parameters and is assumed to be continuous, monotonically increasing, 
bounded by a constant amax, and to satisfy $(s) < s (see Figure 1). 

If M is an upper bound for ||F|^](x)||, then Vl|F[fc] ^ 1/^^ 
(cf. Lemma 3.2). Hence the relaxation function $rei needs only be defined 
on [l/M^,oo). In particular, if one chooses <I>i.ci(s) = amin being constant 
on that interval, then dk ~ ^niin 

and the L-SDK method reduces to the 
loping Landweber-Kaczmarz (l-LK) method considered in [10, 11]. The 



2 




Figure 1: Typical examples for relaxation function $rei- 

convergence analysis of the L-LK method requires amin ^ 1/M^, whereas 
the adaptive choice of the relaxation parameters in the present paper allows 
ak being much larger than 

The L-SDK method consists in incorporating the Kaczmarz strategy 
(with the loping parameters uJk) in the steepest descent method. This 
strategy is analog to the one introduced in [11] regarding the Landweber- 
Kaczmarz iteration. As usual in Kaczmarz-type algorithms, a group of N 
subsequent steps (starting at some multiple k of N) shall be called a cycle. 
The iteration should be terminated when, for the first time, all are equal 
within a cycle. That is, we stop the iteration at 

fcf := argmin{/iV G N : xf^ = xf^+i = • • • = xf^+^.J , (8) 

Notice that /cf is the smallest multiple of N such that 

= = • • • = Xl^s_^_]\f_l . (9) 

In the case of noise free data, (^j = in (1), we choose Wfc = 1 and the iteration 
(4) - (7) reduces to the Steepest-Descent-Kaczmarz (SDK) method, which is 
closely related to the Landweber-Kaczmarz (LK) method considered in [17]. 

It is worth noticing that, for noisy data, the L-SDK method is fundamen- 
tally different from the SDK method: The bang-bang relaxation parameter 
UJk effects that the iterates defined in (4) become stationary if all components 
of the residual vector ||Fj(a;^) —yf\\ fall below a pre-specified threshold. This 
characteristic renders (4) - (7) a regularization method (see Section 3). An- 
other consequence of using these relaxation parameters is the fact that, after 
a large number of iterations, tuk will vanish for some k within each iteration 
cycle. Therefore, the computational expensive evaluation of Fj^j(xfc)* might 
be loped, making the L-SDK method in (4) - (7) a fast alternative to the LK 
method in [17]. Since in praxis the steepest descent method performs better 
than the Landweber method, the L-SDK is expected to be more efficient 
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than the L-LK method [10, 11]. Our numerical experiments (mainly for the 
nonhnear problem considered in Section 5) corroborate this conjecture. 

The article is outlined as follows. In Section 2 we formulate basic as- 
sumptions and derive some auxiliary estimates required for the analysis. In 
Section 3 we provide a convergence analysis for the L-SDK method. In Sec- 
tions 4 and 5 we compare the numerical performance of the L-SDK method 
with other standard methods for inverse problems in photoacoustic tomog- 
raphy and in semiconductors respectively. 



2 Assumptions and Basic Results 

We begin this section by introducing some assumptions, that are necessary 
for the convergence analysis presented in the next section. These assump- 
tions derive from the classical assumptions used in the analysis of iterative 
regularization methods [7, 16, 29]. 

First, we assume that the operators Fj are continuously Frechet differ- 
entiable, and also that there exist xq & X , M > 0, and p > such that 

N-l 

||F:(x)|| <M, xeBpixo) C n A. (10) 

i=0 

Notice that Xq = xq is used as starting value of the L-SDK iteration. Next 
we make an uniform assumption on the nonlinearity of the operators Fj. 
Namely, we assume that the local tangential cone condition [7, 16] 

\\F,{x)-F,{x)-F[{x){x-x)\\y ^^^^ 

< T/||Fj(x) - Fi(x)||y , x,x£Bp{xo) 

holds for some rj < 1/2. Moreover, we assume the existence of and element 

X* G Bp/2{xQ) such that F(x*) = y . (12) 

where y = {yo, . . . , yN-i) are the exact data satisfying (1). 

We are now in position to choose the positive constants amin) t in (7), 
(6). For the rest of this article we shall assume 

amin := ^rcl (l/M^) , (13) 

r > > 2. (14) 

- l-2r] - ^ ' 

In particular, for linear problems we can choose r equal to 2. 

In the sequel we verify some basic results that are necessary for the 
convergence analysis derived in the next section. The first result concerns 
the well-definedness and positivity of the relaxation parameter a^- 
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Lemma 2.1. Let assumptions (10) - (12) be satisfied. Then the coefficients 
Uk in (7) are well-defined and positive. 

Proof. If ujk = 0, the assertion follows from (7). If ujk = 1, then ||F[fc](x|) — 
yfk]\\ ^ T^[k] ^-iid the assertion is a consequence of [29, Lemma 3.1], applied 
to F[^.] instead of F. □ 

In the next lemma we prove an estimate for the step size of the L-SDK 
iteration. 

Lemma 2.2. Let Sk and at he defined by (5) and (7). Then 

akWskf < ||F[fc](xi)-yf,]f, keN. (15) 
Proof. It is enough to consider the case Wfc = 1. It follows from (7) that 

= (piU^) - PiSbp ■ 

Moreover, from the definition of we obtain 

\\F[k]{xi)sk\\ = ||F(,](4)Fjfc](4)*[F[fc](xf)-4]]||, 

\\skf < ||F{,](xi)F(,](4)*[F[,](4) - ||F[,](4) - 4,11 . 

Now, substituting the last two expressions in (16), shows (15). □ 

The following Lemma is an important auxiliary result, which will be used 
at several places throughout this article. 



Lemma 2.3. Let xf., ak, ojk, and Sk be defined by (4) - (7) and assume that 
(10) - (12) hold true. Lf £ Bp/2{x*) for some k > 0, then 

114+1 -x*f-||4-x*f (17) 

< a;fcafc||F[,](4) - yf,]!! ((2ry - l)||F[fc](xf ) - yf.^W + 2(1 + 77)5[,]) . 

Proof. If iOk = 0, then x^+i = x^ and (17) follows with equality. If = 1, 
it follows from (4) and (5) and Lemma 2.2 that 

ll-^fc+l X II ||X^ X II 

= - X*, x^+i - xf ) + ||x^_,_i - xf,|p 
= 2ak{xi - X*, F|,](xi)*(4] - + 

< 2tofc(4] - F[fc](xiO,F|,](xt)(xt - X*)) + afc||F[fc](xt) - yf.jf 

< afc(2(4] - F[,](xf ), F|,](x^)(xi - x*) - F[,](x*) + F[,](xiO) 

+ 2(4.] - F[fc](xi), y[,] - - - F[,](xf )f ) . 
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Now, applying (11) with x = x* and x = G Bpi2{x*) C Bp{xo), leads to 

„*l|2 ||^<5 „*||2 



< a;fcafc||F[fc](4) - yf,]|| ^277||F[fc] (x^^) - + 25[k] - \\F[k]ixi) - yf^jl 
The last inequality and (1) show (17). □ 

Our next goal is to prove a monotony property, known to be satisfied 
by other iterative regularization methods, e.g., by the Landweber [7], the 
steepest descent [29], the LK [17], and the L-LK [11] method. 

Proposition 2.4 (Monotonicity). Under the assumptions of Lemma 2.3, 

\\xi_^i- x*f <\\xi- x*f , ken. (18) 

Moreover, all iterates xf, remain in Bp/2{x*) C Bp{xo) and satisfy (17). 

Proof. From (12) it follows that xq G Bpj2{x*). If uj^ = 0, then xi satisfies 
(18) with equality and xi G Bp/2{x*) C Bp{xo). If co^ ^ 0, then Lemma 2.3 
implies 



|x^ - x*||2 > (27? - l)llFo(xg) - /'°|| + 2(1 + r?)50 
> (^of(2r?- l)r + 2(l + r/)) > 0. 



(5 *||2 
Xi — X II 



Therefore (18), for A; = 0, follows from (14). In particular, x\ G Bp/2{x*). 
An inductive argument implies (18) and that Xk G Bpi2{x*) C i3p(xo) for 
all /c G N. The assertions therefore follows from Lemma 2.3. □ 



3 Convergence Analysis of the Loping Steepest 
Descent Kaczmarz Method 

In this section we provide a complete convergence analysis for the L-SDK 
iteration, showing that it is a convergent regularization method in the sense 
of [7] (see Theorems 3.3 and 3.6 below). Throughout this section, we assume 
that (10) - (14) hold, and that x^, a^, w^, and Sk are defined by (4) - (7). 

Our first goal is to prove convergence of the L-SDK iteration for (5 = 0. 
For exact data y = {yo, . . . ,yj\[~i), the iterates in (4) are denoted by x^.^ 

Lemma 3.1. There exists an XQ-minimal norm solution of (2) in Bp/2{xo), 
i.e., a solution x^ of (2) such that 

\\x^ — xqW = inf { ||x — xoll : X G -Bp/2(a^o) ^^^'d F(x) = y} . 

Moreover, x^ is the only solution of (2) in Bp/2{xo) H (xq + ker(F'(x^))"'-) . 

^This is a standard notation used in the literature. 
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Proof. Lemma 3.1 is a consequence of [13, Proposition 2.1]. A detailed proof 
can be found in [16]. □ 

Lemma 3.2. For all € N, we have Uk > amin- 

Proof. For = the claimed estimate holds with equality. If Wfc = l, it 
follows from (10) that 

PfciiViiFj,,](4)«fcf > w[kH)\\-^ > m'- 

Now the monotonicity of $rei implies ak > <I>rei(-^~^) = amin- D 

Throughout the rest of this article, denotes the XQ-minimal norm 
solution of (2). We define := — x^. Prom Proposition 2.4 it follows 
that (17) holds for all k. By summing over all k, this leads to 

°° II tl|2 

J]a.||y[,]-F[,](rr,)f < "^""^ " < oo . (19) 

Equation (19) and the monotony of ||efc|| shown in Proposition 2.4 are main 
ingredients in the following proof of the convergence of the SDK iteration. 

Theorem 3.3 (Convergence for Exact Data). For exact data, the iteration 
(xfc) converges to a solution of (2), as k ^ oo. Moreover, if 

M{F'{x'^)) C J\f{F'{x)) for all x G Bp{xo) , (20) 

then — > x^. 

Proof. From (18) it follows that ||efc|| decreases monotonically and therefore 
that ||efc|| converges to some e > 0. In the following we show that is in 
fact a Cauchy sequence. 

For k = koN + ki and I = loN + h with k < I and ki,li € {0, . . . , N -1}, 
let uq £ {ko, . . . , /q} be such that 

N-l N-1 



\\Fij^{xNno+ii) - VhW < \\Fii{xNio+ii) - VhW , io e {^0, • • • ,^o} • 

(21) 



n=0 n=0 



Then, with n := Nuq + N — I, we have 

Ijcfc - ei\\ < \\ek - en\\ + \\ei - e„|| (22) 

and 



I _ ||2 _ II ||2 _ II ||2 _i_ 9/ _ 
l^^n — ll^fcll I Infill I ^\^n Cfcjfi. 

II I|2 II ||2 II ||2 I o/ 

||en-ei|| = \\ei\\ - \\en\\ +2(en-ez,e 



(23) 



n/ ) 
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For k — > oo, the first two terms of (23) converge to e — e = 0. Tlierefore, 
in order to sliow that is a Cauchy sequence, it is sufficient to prove that 
(cn — Cfc, Bn) and {cn — €[, Bn) Converge to zero as A; ^ oo. 

To that end, we write i = Niq + ii, ii £ {0, . . . , N — 1} and set i* := 
NriQ + ii. Then, using the definition of the steepest descent Kaczmarz 
iteration it follows that 

I (en - efc,en)| 

n-1 



i=k 
n-1 

- \{yi^ ~ Fn(a;i),F-^(xi)(x''' - Xi*) + ¥[^{xi){xi* - Xn)) 

i=k 
n-1 



n— 1 



+ - Fn(x^)|| ||F^,(xi)(xi. - x„)|| (24) 

i=k 



From (11) it follows that 

||Faxi)(xt-x,*)l| <2(l+r?)||2/i, -Fi,(xi)|| + (l + r?)||yii -Fi,(xi.; 



(25) 



Again using the definition of the steepest descent Kaczmarz iteration and 
equations (7), (10), it follows that 

||F-^(Xi)(Xi*-X„)|| < M ||xj. - Xn\\ 
N-2 

< M E aj\\F'j{xNno+j)* {Fj{xNno+j) - yj)\\ 

j=ii 

N-1 

< amaxM^ \\Fj{xNno+j) " VjW ■ (26) 

3=0 

Substituting (25), (26) in (24) leads to 

I (Cn Cfc , Gji) I 

n-1 N-1 /N-1 

< c E E WVi, -Fi^{xmo+ii)\\ Yl \\^i(^Nno+j) - VjW 
io=koh=0 \j=0 
n-1 N-1 



< C E (E llyii - Fi^{xNio+h)\\) 
io=ko ii=0 
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with c := amax(3 + 3r] + ckmax-^^)- Here we made use of (21). So, we finally 
obtain the estimate 

Nc ^"^ 

|(e„ - efc,e„)| < ^ amo+iAlyh - Fi^{xmo+ii)f ■ 

iO=ko ii=U 

Because of (19), the last sum tends to zero for k = Nko + ki — > oo, and 
therefore {en,en — e^) 0. Analogously one shows that {en, en — e/) — > 0. 
Therefore is a Cauchy sequence and = — converges to an element 
X* G X. Because all residuals ||F[fc](3;fc) — tend to zero, x* is solution 
of (2). 

Now assume Af(F'{x^)) C A^(F(x)), for x G Bp^xo). Then from the 
definition of it follows that 

Xk+i -Xke 7^(F^,,](xfc)*) C AA(Fj;-](rEfc))^ C M{F'{xk))^ C Af{F'{x^))^ . 

An inductive argument shows that all iterates Xk are elements of xq + 
AA(F'(x^ ))"*-. Together with the continuity of F'(x^^) this implies that x* E 
a;o+AA(F'(xt))-L. By Lemma 3.1, x^ is the only solution of (2) in Bp/2{xo)^ 
[xo + Af(F'{x^))-^), and so the second assertion follows. □ 

The second goal in this section is to prove that x^^ converges to a solution 
of (2), as 5 ^ 0. First we verify that, for noisy data, the stopping index k^ 
defined in (8) is finite. 

Proposition 3.4 (Stopping Index). Assume 5mm := min{(5o, . . . , ^at-i} > 
0. Then kl defined in (8) is finite, and 

\\F^{xis)-y!\\<r6i, i = 0,...,N-l. (27) 

Proof. Assume that for every I G N, there exists i{l) G {0, . . . , — 1} such 
that x;7v+i{i) 7^ xiN. From Proposition 2.4 follows that we can apply (17) 
recursively for k = 1, . . . , IN and obtain 

IN 



-\\xo-x*f <^u;fcafc||F[fc](xf,) 

k=l 

(2(l + r?)5[,]-(l-2r/)||F[,](xi)-4,]||) , ZgN. 
Using the fact that either a;^ = or ||F[^.](x^) — > rd^k], we obtain 

IN 

\\xo-x*f > (r(l - 2r/) - 2(1 + 7?)) ^ uJkak6^,^\\F^kM) - yfk]\\ ■ (28) 



k=l 
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Equation (28), Lemma 3.2 and the fact that xii / xiin for all /' G N, 
imply 

lko-x*f > (r(l-2r?)-2(l + ?7))/a^in5,ni„(r5n,in), I en. (29) 

The right hand side of (29) tends to infinity, which gives a contradiction. 
Consequently, G N : xi^j^i = xin ,0<i<N — 1}^^ and the infimum 
in (8) takes a finite value. 

To prove (27), assume to the contrary, that ||Fj(x^a) — yf\\ > rdi for 
some i G {0, . . . , — 1}. From (6) and (8) it follows that, ijJj^s = 1 and 
— ^^ki+i+i respectively. Thus, Proposition 2.4 and Lemma 2.1 imply 

< (2r/ - l)||F,(x^,) - II + 2(1 + r])5i < Si ((2r? - 1)t + 2(1 + r?)) . 

This contradicts (14), concluding the proof of (27). □ 

The last auxiliary result concerns the continuity of at 5 = 0. For 
y, G Y^, 5 > 0, and A; G N we define 

Ak{5,y,y') := ^fcFj,](4)*(F[fc](xt) - yf,,) - F|,](xfc)*(F[fc](xfe) - y[fc]) . 

Lemma 3.5. For all A; G N, 

lim sup|||Afc(5,y,y^)|| : G Y^, \\yi - yf\\ < = . (30) 

Moreover, xl._^-^ — > x/^^i, as (5 — > 0. 

Proof. We prove Lemma 3.5 by induction. The case A; = is similar to the 
general case and is omitted. 

Now, assume k > and that (30) holds for all k' < k. First we note that 
(30) and the continuity of <^rei obviously imply xf.^^ — > x^+i, as 5^0. For 
the proof of (30) we consider two cases. In the first case, = 1, we have 

\\Ak{5,y,y')\\ = ||F{,](xi)*(F[,](xf ) - yf,]) - F|,](xfc)*(F[,](xfc) - y[fc])|| . 
In the second case, Wfc = 0, we have ||F[fc] (xf.) — y^^j || < t6'' and consequently 

\\A,i5,y,y')\\ < ||Fj,](x,)*(F[,](xfc) - y[,,])|| 

< l|F|fc](4)ll (||F[,,](xfc) -F[fc](xi)|| + ||F[fc](4) -yfll + U - y[k]\\) 

< llFjfc] (4)11 ( l|F[fc] (xfc) - F[,] (4)11 + (r + 1) 6[k]) ■ 

Now (30) follows from (10), the continuity of -Fj^] and and the induction 
hypothesis (which implies xf. — > Xk 
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Theorem 3.6 (Convergence for Noisy Data). Assume {5q, . . . ,5j^_^) is a 

sequence in (0, oo)^ with limj^oo^l = 0- Let {uq, . . . ,y-'^_^) he a sequence 
of noisy data satisfying 

H-yi\\<4^ i = 0,...,iV-l, jGN, 

and let := k^:{5^,y^) denote the corresponding stopping index defined in 
(8). Then Xy converges to a solution of (2), as j ^ oo. Moreover, if (20) 
holds, then — > x^. 

Proof. Let x* denote the limit of the iterates Xk which is a solution of (2), 
cf. Theorem 3.3. Prom Lemma 3.5 and the continuity of Fj we know that, 
for any fixed k gN, 

xf^Xk, Fi{xf)^Fi{xk), asj^oo. (31) 

To show that x^^j — > x*, we first assume that k^ has a finite accumulation 
point fc^,. Without loss of generality we may assume that k^ = k^ for all 
j G N. From Proposition 3.4 we know that \\yf — Fj(x^'' )|| < rSj and, by 
taking the limit j — > oo, that = Pj(xfcJ. Consequently Xk^ = x* and 

It remains to consider the case where A:-' — > oo as j — > oo. To that end 
let e > 0. Without loss of generality we assume that k^ is monotonically 
increasing. According to Theorem 3.3 we can choose n G N such that 
llxfcn — < e/2. Equation (31) implies that there exists jo > n such 
that \\xf„ — Xkn\\ < e/2 for all j > jq. This and Proposition 2.4 imply 

W-^y ~ ^ II — iFfc" ~ ^ II 

< llxfn - Xfcn|| + ||xfcn - < | + | = £, for j > Jq • 

Consequently, Xy — > x* . 

If (20) holds true, then by Theorem 3.3, x* = x^. Therefore x^ x\ 
which concludes the proof. □ 

Remark 3.7. In standard iterative regularization methods the number of 
performed iterations plays the role of the regularization parameter [7, 16]. A 
parameter choice rule corresponds to the choice of an appropriate stopping 
index fef = k{5,y^). 

For the loping Kaczmarz iterations analyzed in this article, the situa- 
tion is quite different. If k is fixed, then the iterates xf., do not depend 
continuously on data yf. However, for a fixed sequence {uJk) of loping pa- 
rameters, the iterates xf. do depend continuously on yf: Now, the loping 
sequences (uJk) play the role of the regularization parameters and the partic- 
ular sequence uJk = ^k{^-,y^), depending on 6i and the noisy data yf, is the 
a-posteriori parameter choice rule. 
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4 Limited View Problem in Photoacoustic Com- 
puted Tomography 

In this section we compare the numerical performance of loping Kaczmarz 
methods applied to a system of linear equations related to a limited view 
problem in photoacoustic computed tomography [8, 18, 28, 34]. 

Let X := L^{D) denote the Hilbert space of all square integrable func- 
tions in the unit disc D C M?, and let Y denote the Hilbert space of all 
functions y : [0,2] ^ M with := y{t)tdt < oo. We consider the 

system 

MiX = yi, i = 0,...,N-l, (32) 

where Mi : X ^Y, 

{Mix)it) := ^ [ x{Ci + ta) dn{a) , t £ [0,2] , (33) 

correspond to a scaled version of the circular mean Radon transform. Solv- 
ing (32) is the crucial step in three-dimensional photoacoustic computed 
tomography with integrating linear detectors [9, 28], where the centers of 
integration, ^j, correspond to the positions of the linear detectors. We are 
particularly interested in the incomplete data case (limited view problem), 
where the centers = [sin{iri/[N — 1)) cos{iri/{N — 1))) are uniformly dis- 
tributed on the semicircle S\ :={£, = , S,'^) £ dD : ^} > 0}. Micro- local 
analysis predicts, that if the centers do not cover the whole circle, certain 
details (the invisible boundaries) of x outside the detection region (convex 
hull of S\) cannot be recovered [21, 27, 35]. 

The operators Mj are linear, bounded, and satisfy || Mj || < 1, [10]. 
For linear bounded operators, the tangential cone condition (11) is satisfied 
with r/ = 0. Consequently, the analysis of Section 3 applies, and the L-SDK 
method (4) - (7) provides a convergent regularization method for solving 
(32). The adjoint of Mj, required in (5) is given by (M*y)(^) = y{\£,i — 
[10]. 

In the following numerical examples, we consider the L-SDK method 
with either the choice $rci(s) = ^^^(aminS, 2) or <^rei(s) = amin (which 
corresponds to the L-LK method). In both cases we use amin = 0.4 or 
c^min = 1) and assume = 50 measurements. The phantom x^, shown 
in the left picture in Figure 2, consists of a superposition of characteristic 
functions and one Gaussian kernel. Data yi = Mj were calculated via 
numerical integration with the trapezoidal rule and 4% noise was added, 
such that — yf||/||yj|| ~ 0.04. The the regularized solutions x^^ with 
amin = 0.4 are depicted in Figure 3. For both, the L-SDK L-LK method, 
all visible parts of the phantom are reconstructed reliable. 

Figure 4 and Figure 5 show the number of actually performed iterations 
and the reconstruction error := ||x| — respectively. For comparison 
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Figure 2: The left picture shows the phantom x^, where the white dots 
indicate the locations of the detectors. The corresponding data {yf)i are 
depicted on the right. 



purposes, the error for the SDK and the LK iteration (without loping pa- 
rameter) are also included. In all cases, the smaller relaxation parameter 
gives the smaller reconstruction errors. This behavior is typically for 



a 



the application of Kaczmarz type iterations to Radon transforms [5, 25]; 
therefore in praxis often relatively small relaxation parameters are chosen. 
For amin = Ij the loping strategy significantly reduces the reconstruction 
error of the no-loping iterations. Also, for amin = 0.4, the regularized so- 
lution of the loping Kaczmarz methods (automatically stopped according 
to (8)) have errors comparable to the optimal solution of their non- loping 
counterparts when stopped after the cycle with minimal error (which is not 
available in practice). 





Figure 3: Numerical reconstructions with a^, 
depicted in Figure 2. 



0.4 of the phantom 
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Figure 4: The x-axis shows the number of cycles, while the number of 
actually performed iterations within each cycle is shown at the y-axis. 
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Figure 5: Evolution of the relative error In 



To point out the effectiveness of the loping Kaczmarz methods for solv- 
ing linear inconsistent systems we included the reconstruction error for the 
CGNE iteration (conjugate gradient [14, 31] applied to normal equations). 
If stooped appropriately the CGNE method is known to be a regularization 
method [7, 12]. As can be seen in Figure 5 the reconstruction error for the 
L-SDK and the L-LK methods is much smaller that that for the CGNE it- 
eration. In Table 1 run times for reconstructing an image on a 120 x 120 grid 





Cycles 


Runtime (sec) 


Error (%) 


L-SDK 


5 


21.9 


18.2 


L-LK 


6 


21.4 


18.5 


SDK 


4 


24.5 


18.2 


LK 


5 


16.9 


18.1 


CGNE 


5 


38.2 


21.6 



Table 1: Comparison of the performance of different iterative methods. The 
non-loping iterations are stopped after the cycle with minimal error. 
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are compared (with non-optimized Matlab implementation on iMac with 2 
GHz Intel Core Duo processor). 

5 An Inverse Doping Problem 

In this section we present another comparison of the numerical performance 
of the L-SDK, L-LK and LK methods. This time we consider an application 
related to inverse doping problems for semiconductors [2, 10, 20, 3] For 
details on the mathematical modeling of this inverse problem we refer the 
reader to [10, Section 3]. 

In what follows we describe the abstract formulation in Hilbert spaces 
of the problem (the so called inverse doping problem in the linearized unipo- 
lar model for current flow measurements). Let O := (0, 1) x (0, 1) C 
be the domain representing the semiconductor device (a diode). The two 
semiconductor contacts are represented by the boundary parts: 

To := {(s, 0) : s e (0, 1)} , Ti := {{s, 1) : s G (0, 1)} , 

(we denote OQd '■= ToUTi) while the insulated surfaces of the semiconductor 
are represented by d^jy := {{0,t) : t E (0,1)} U {(l,t) : t E (0,1)}. 
This specific inverse doping problem can be reduced to the identification of 
the positive parameter function x (the doping profile C is related to x by 
C = X — A^A(lnx)) in the model 

A^n V • (x(0 Vn) = , in (34) 
u = U{0, ondQD (35) 
Vn-zv = 0, on d^N (36) 

from measurements of the Voltage-Current map (the forward operator) 

which maps an applied potential U at d^l^ to the corresponding total current 
flow T,x{U) through the contact Ti. Here fin, A are positive constants and 
Vbi is a known logarithmic function defined on dQ^- 

Due to the nature of the practical experiments that can be performed 
on a factory environment, some restrictions on the data have to be taken 
into account: 

1. The voltage profiles U E H'^/'^idVLo) must satisfy U{i) = at the 
contact Ti. 
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Figure 6: In the top left pieture, the doping profile to be identified. In the top 
right picture, a typical voltage profile Ui and the corresponding solution u of (34) 
- (36). The initial guess used for the L-SDK, L-LK and LK iterative methods is 
shown in the bottom picture. The boundary parts Tq and Fi correspond to the top 
right and to the lower left edge respectively (the origin is the right corner). 

2. The parameter x has to be determined from a finite number of mea- 
surements, i.e. from the data 

:=S,(^7,) ey :=M, i = 0,...,N-l, (37) 

where the Ui G H^^'^{dO,£)) are prescribed voltage profiles satisfying 
Item 1. 

Therefore we can model the inverse doping problem with a system of oper- 
ator equations of the form (2), namely 

Fi{x)=y!, i = 0,...,N-l, 

where x G L^(J7) =: X is the unknown parameter, G M. =: Y denote the 
measured data, Fi : X ^ Y defined by Fi{x) := Sa;(C/j) are the parameter 
to output maps, with domains of definition 

Di := {x G L°°(r2) : < Xmin < x < Xmax, a.e.} . 

It is worth mentioning that, although the operators Fj are Frechet differen- 
tiable, they do not satisfy the tangential cone condition (11). Therefore, the 
convergence results derived in Section 3 cannot be applied. 
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Figure 7: Comparison between the L-SDK, L-LK and LK mctliods. The top two 
pictures show the iterative errors obtained by the L-SDK iteration after 10 and 81 
cycles. The two central pictures show the iterative errors obtained by the L-LK 
iteration after 80 and 121 cycles. The two bottom pictures show the iterative errors 
obtained by the LK after 120 and 380 cycles. 



In the following numerical examples we assume that = 11 Dirichlet- 
Neumann pairs {Ui,Fi{x')) of measurement data are available. The fixed 
inputs Ui, are chosen to be piecewise constant functions supported in Tq, 

TT t \ / 1 \s- Si\<h 

■= [ else 

where the points Sj are uniformly distributed on Tq and h = 1/32. The 
doping profile to be reconstructed is shown in Figure 6 (top left picture). 
The top right picture of Figure 6 shows a typical voltage profile Uj (applied 
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at To) as well as the corresponding solution u of (34) - (36). In these pictures, 
as well as in the forthcoming ones, Fi appears on the lower left edge and Tq 
on the top right edge (the origin corresponds to the upper right corner). 

In Figure 7 we show the evolution of the iteration error for the L-SDK, 
L-LK and LK methods. The same initial guess was used for the three 
methods (see Figure 6). In our computations we chose r = 2.5 in (14). The 
stopping rule for the L-SDK method is satisfied after 81 cycles. For the 
L-LK method, the same stopping criteria is reached only after 121 cycles. 
In order to obtain the same accuracy with the LK method, 380 cycles are 
required. In the top pictures of Figure 7 one can see the iteration error for 
the L-SDK method after 10 and 81 cycles. For comparison purposes, the 
iteration error for the L-LK method is shown after 80 and 121 cycles (see 
the central pictures of Figure 7). The bottom pictures of Figure 7 show the 
iteration error for the LK method after 120 and 380 cycles. The number of 
actually computed iterative steps within each cycle of the L-SDK and L-LK 
methods is shown in Figure 8. 

As one can see in Figure 8, no more than 2 steepest descent steps per 
cycle are computed after the 14-th cycle of the L-SDK method. Analogously, 
no more than 2 Landweber steps per cycle are computed after the 37-th cycle 
of the L-LK method. In total, for the computation of the LK-approximation 
in Figure 7 (380 cycles), 4180 Landweber steps are needed, while the L- 
LK-approximation (121 cycles) requires the computation of 258 Landweber 
steps and the L-SDK- approximation (81 cycles) requires the computation 
of 184 steepest descent steps. The L-LK method requires almost 50% more 
cycles than the L-SDK method in order to reach the stopping criteria (8). 
Moreover, the LK method requires almost three times more cycles than 
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Figure 8: Comparison between the performance of L-SDK and L-LK methods. 
The solid (blue) line shows the actually performed number of steps within each 
cycle of the L-SDK method, while the dashed (red) line gives the corresponding 
information with respect to the L-LK method. 
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the L-LK method in order to achieve the same accuracy (see [10] for other 
comparisons between the LK and L-LK methods). 

The efficiency of the L-SDK method becomes even more evident when 
we compare the total number of actually performed iterative steps. Each 
cycle of the LK method requires the computation of 11 steps, while in the 
L-SDK and L-LK methods the number of actually performed steps per cycle 
is very small after a few number of cycles. 

6 Conclusions 

In this paper we propose a new iterative method for inverse problems of the 
form (2), namely the L-SDK method. As a by-product we also formulated 
the SDK iteration, which is the steepest descent counterpart of the LK 
method [17]. In the L-SDK iteration we omit an update of the SDK iteration 
(within one cycle) if corresponding i-th residual is below some threshold. 
Consequently, the L-SDK method is not stopped until all residuals are below 
the specified threshold. We provided a complete convergence analysis for the 
L-SDK iteration, proving that it is a convergent regularization method in 
the sense of [7] . 

The abstract theory was applied to thermoacoustic computed tomogra- 
phy and an inverse problem for semiconductors. In both applications the 
L-SDK method turned out to be an efficient iterative regularization method. 
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